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1 Program Summary 

Title of program: 
Catalog number: 
Program obtainable 

Computer on which the program 
has been thoroughly tested: 

Operating system: 

Programming language used: 

Memory required to execute 

No. of bits in a word 

Has the code been vectorized? 

Number of lines 

in distributed program: 

Other programs used in SOPHIA 
in modified form 



Nature of physical problem: 



SOPHIA 2.0 
from A. Miicke 

e-mail: amuecke@physics.adelaide.edu.au 

DEC- Alpha and Intel-Pentium based 
workstations 

UNIX, Linux, Open- VMS 

FORTRAN 77 

<1 megabyte 

64 

no 

16500 



Rndm 

Processor independent random number 

generator based on Ref. [0] 

Jetset 7.4 

Lund Monte Carlo for 

jet fragmentation 

DECSIB 

SiBYLL routine which decays 
unstable particles 

Simulation of minimum bias photo- 
production for astrophysical applications 
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Method of solution: Monte Carlo simulation of individual events 

for given nucleon and photon energies, 
photon energies are sampled 
from various distributions. 

Incident ambient photon fields limited 
to power law and blackbody spectra so far. 
No tests were done for center-of-mass 
energies ^/s > 1000 GeV. 

Typical running time: 10000 events at a center-of-mass energy 

of 1.5 GeV require a typical 
CPU time of about 75 seconds. 



Restrictions on the complexity 
of the problem 



3 



2 Introduction 



The cosmic ray spectrum extends to extremely high energies. Giant air showers have been 
observed with energy exceeding ~ 10^^ GeV 0. Energy losses due to interactions with 
ambient photons can become important, even dominant for such energetic nucleons, above 
the threshold for pion production. Photoproduction of hadrons is expected to cause a dis- 
tortion of the ultra-high energy cosmic ray (CR) spectrum by interactions of the nucleons 
with the microwave background (the Greisen-Zatsepin-Kuzmin cutoff P, ^; see also [§ for 
additional references), but it may also be relevant to the observed high energy gamma ray 
emission from jets of Active Galactic Nuclei (AGN) (e.g. 0, [1^) or Gamma- Ray Bursts 
(GRB) [|TT|. Moreover, it is the major source process for the predicted fluxes of very high 



energy cosmic neutrinos (e.g. |T2|, |T3|, 0, |T^, [16[|). 



The photohadronic cross section at low interaction energies is dominated by the A(1232) 
resonance. Since the low energy region of the cross section is emphasized in many astrophys- 
ical applications, the cross section and decay properties of the prominent A-resonance have 
often been used as an approximation for photopion production, and the subsequent produc- 



tion of gamma rays and neutrinos ||T^, 0. As discussed in [|1^, this approximation is 
only valid for a restricted number of cases, and does not describe sufficiently well the whole 
energy range of photohadronic interactions. A more sophisticated photoproduction simula- 
tion code is needed to cover the center-of-mass energy range of about y/s ^ 1 — 10'^ GeV, 
which is important in many astrophysical applications. 

In this paper we present the newly developed Monte-Carlo event generator SOPHIA 
(Simulations Of Photo Hadronic Interactions in Astrophysics), which we wrote as a tool 
for solving problems connected to photohadronic processes in astrophysical environments. 
The philosophy of the development of SOPHIA has been to implement well established phe- 
nomenological models, symmetries of hadronic interactions in a way that describes correctly 
the available exclusive and inclusive photohadronic cross section data obtained at fixed tar- 
get and collider experiments. 

The paper is organized as follows. After introducing the kinematics of A^7-interactions 
(Sect. 3) we give a brief physical description of the relevant photohadronic interaction pro- 
cesses (Sect. 4). The implementation of these processes into the SOPHIA event generator, 
together with the method of cross section decomposition and parametrization, is described 
in Sect. 5. The structure of the program is outlined in Sec. 6, and a comparison of the model 
results with experimental data is provided in Sec. 7. The definitions of special functions 
and tables of parameters used in the cross section and final state parametrization, as well 
as a compilation of the important routines and functions used in the code, are given in the 
appendices. 

Unless noted otherwise, natural units {h = c = e = 1) are used throughout this paper, 
with GeV as the general unit. In this notation, cross sections will be in GeV~^. A general 
exception is Section 4, where numerical parametrizations of cross sections are given in /xbarn. 
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The relevant conversion constant is {Key = 389.37966 GeV-Vbarn. 

3 The physics of photohadronic interactions 
3.1 Kinematics of collisions 

There are three reference frames involved in the description of an astrophysical photohadronic 
interaction: (i) the astrophysical lab frame (LF), (ii) the rest frame of the nucleonQ (NRF), 
and (iii) the center-of-mass frame (CMF) of the interaction. For example, in the LF, the 
initial state can be characterized by the nucleon energy En, the photon energy e, and the 
interaction angle 6 

cos9 = {pn ■ p^)/pNENe ■ (1) 

where and p^ denote the nucleon and photon momenta. The Lorentz factor is = 
En/uin = (1 — Pn)"^^'^ with rriN being the nucleon mass. The corresponding quantities in 
the NRF and CMF are marked with a prime (') and an asterisk (*), respectively. Fixed target 
accelerator experiments where a photon beam interacts with a proton target are performed 
in the NRF. In astronomical applications we assume that the LF can be chosen such that 
a the photon distribution function is isotropic. The LF may therefore be different from the 
astronomical observer's frame if, for example the emission region is moving relative to the 
observer such as in AGN jets or GRB. The interaction rate of the nucleon in the LF is given 
by 

1 POO ni^] /"^max 

R{En) = —;rir / * ^ / - ^1>N^{S) , (2) 

where a^-y is the total photohadronic cross section and 

s = m\ + 2i?Are(l — cos 6) = m\ + 2mj\[e', (3) 

is the square of the center-of-mass energy. The lowest threshold energy for photomeson 
production is y^s^ = rriN + 171^^0. The remaining quantities in Eq. (0) are eth = {sth — 
m%)/2{EN + Pn) and s^ax = m]^ + 2Eive(l + I3n)- 

The final state of the interaction is described by a number of possible channels, each 
of which has Nf^c particles in the final state. At threshold, the phase space volume vanishes, 
which requires kinematically for the partial cross section cxc — for s Sth,c = [Si^i]^- 
Above threshold, each final state channel has 3A^/ c — 4 degrees of freedom given by the 3- 
momentum components {PiiXii'Pi) of the particles, constrained by energy and momentum 
conservation. Here p,, x-ii and 0j are the particle momentum, and it's polar and azimuthal 
angles with respect to the initial nucleon momentum, respectively. One of the 0, angles 

^In the following the subscript TV is used if no distinction between protons and neutrons is made. Inter- 
actions for both protons and neutrons are implemented in SOPHIA. 
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can be chosen to be distributed isotropically since we consider only the scattering of unpo- 
larizcd photons and nucleons, all other variables are determined by the interaction physics 
through the differential cross sections. A distinguished role in the final state is played by the 
"leading-baryon" , which is considered to carry the baryonic quantum numbers of the incom- 
ing nucleon. For this particle, the Lorentz invariant 4-momentum transfer t — (Pjv — -Pfinai)^ 
is often used as a final state variable. At small s, many interaction channels can be reduced 
to 2-particle final states, for which da/dt gives a complete description. 



3.2 Interaction processes 

Photon-proton interactions arc dominated by resonance production at low energies. The in- 
coming baryon is excited to a baryonic resonance due to the absorption of the photon. Such 
resonances have very short life times and decay immediately into other hadrons. Conse- 
quently, the Nj cross section exhibits a strong energy dependence with clearly visible reso- 
nance peaks. Another process being important at low energy is the incoherent interaction of 
photons with the virtual structure of the nucleon. This process is called direct meson pro- 
duction. Eventually, at high interaction energies {^/s > 2 GeV) the total interaction cross 
section becomes approximately energy-independent, while the contributions from resonances 
and the direct interaction channels decrease. In this energy range, photon-hadron interac- 
tions are dominated by inelastic multiparticle production (also called multipion production). 



3.2.1 Bciryon resonance excitation and decay 

The energy range from the photopion threshold energy -^^.j^ Ri 1.08 GeV for 7A^-interactions 
up to ^/s ~ 2 GeV is dominated by the process of resonant absorption of a photon by the 
nucleon with the subsequent emission of particles, i.e. the excitation and decay of baryon 
resonances. The cross section for the production of a resonance with angular momentum J 
is given by the Breit-Wigner formula 

where M and F are the nominal mass and the width of the resonance. 6^ is the branching ratio 
for photo-decay of the resonance, which is identical to the probability of photoexcitation. 
The decay of baryon resonances is generally dominated by hadronic channels. The exclusive 
cross sections for the resonant contribution to a hadronic channel with branching ratio be 
can be written as 

(7c(s;M,F, J) = 6,(7bw(s;M,F, J), (5) 

with X^c^c = 1 — ^7 ~ 1. Most decay channels produce two-particle intermediate or final 
states, some of them again involving resonances. For the pion-nucleon decay channel, Ntt, 
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the angular distribution of the final state is given by 



l^oc E IflAdxf , (6) 
ct cos X xZ^j ^ 

where x* denotes the scattering angle in the CMF and // , are the A^vr-helicity amplitudes. 

2 ''^ 

The functions rff ^ire commonly used angular distribution functions which are defined 

A, 2 

on the basis of spherical harmonics. The Nir helicity amplitudes can be determined from 
the helicity amplitudes Ai and Aa for photoexcitation (see Ref. for details), which 



are measured for many baryon resonances [^. The same expression applies to other final 
states involving a nucleon and an isospin-0 meson (e.g., Nrj). For decay channels with other 
spin parameters, however, the situation is more complex, and we assume for simplicity an 
isotropic decay of the resonance. 

Baryon resonances are distinguished by their isospin into A^-resonances (/ = |, as for 
the unexcited nucleon) and A-resonances (/ = |). The charge branching ratios 6iso of the 
resonance decay follow from isospin symmetry. For example, the branching ratios for the 
decay into a two-particle final state involving a A^- or A-baryon and an / = 1 meson (vr or 
p) are given in Table 0. Here A/3 is the difference in the isospin 3-component of the baryon 
between initial and final state (the baryon charge is Q b = h + \) ■ In contrast to the strong 
decay channels, the electromagnetic excitation of the resonance does not conserve isospin. 
Hence, the resonance excitation strengths for py and n'j interactions are not related to each 
other by isospin symmetry and have to be determined experimentally. 



Table 1: Isospin (charge) branching ratios for the decay of a resonance with isospin /res 
and charge /s + | into a final state containing a baryon Bf with isospin /^^ and charge 
+ A/3 + |, and a meson of isospin 1 (tt or p). For example, the decay A^"*" A^+tt^ 
corresponds to /res = |, Isf = -^3 = |) and A/3 = +1, thus /3 A/3 > and 6iso = |- 



^iso 


Ib, = 1/2 {Bf = N) 


Ib, = 3/2 {Bf = A) 




A/3 ^0 A/3 = 


I3 A/3 < A/3 = /3 A/3 > 


1/2 

3/2 


2/3 1/3 
1/3 2/3 


1/6 1/3 1/2 
8/15 1/15 2/5 



3.2.2 Direct pion production 

Direct pion production can be considered as electromagnetic scattering by virtual charged 
mesons, which are the quantum-mechanical representation of the (color-neutral) strong force 
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field around the baryon. The interacting virtual meson gains enough momentum to mate- 
rialize. Experimentally the direct production of charged pions is observed as a relatively 
structureless background in the iVvr^ and Avr''^ final states in photon-nucleon interactions. 

In terms of Feynman graphs, this process is represented by the t-channel exchange of 
a meson. Here, t is the squared 4-momentum transfer from the initial to the final state 
baryon. The graph has a strong vertex at the baryon branch and an electromagnetic vertex 
for the photon interaction. At the strong vertex, the baryon may be excited and change its 
isospin. Isospin combination rules determine the iso-branching ratios in the same way as for 
resonance decay (Table |l| for Jj-es = |). The presence of the electromagnetic vertex requires 
that the particle the photon couples to is charged. Thus direct processes with A/3 = 
branches (e.g. 'jp —>■ n^p) are strongly suppressed. 

The low energy structure of the direct cross section is not well constrained. At high 
energies, Regge theory of the pion exchange implies that crdir(s) oc [|^, |2^. The angular 
distribution of the process is strongly forward peaked and can be parametrized for small |t| 
by ^ 

oc exp(6dirt) • (7) 
with an experimentally determined slope of 6dir ~ 12 GeV~^ 



The total cross section for a direct scattering process is roughly oc , where rrit is the 
(nominal) mass of the exchanged virtual particle. Therefore, the direct production of pions is 
dominant, while the contributions from the exchange of heavier mesons are suppressed. The 
same applies to direct reactions which involve the exchange of a virtual baryon (u-channel 
exchange). However, with increasing energy, more and more channels add to the direct cross 
section, and this makes an explicit treatment difficult. 

3.2.3 High energy processes 

Phenomenologically, high energy interactions can be interpreted as reggeon and pomeron 
exchange processes. Both the reggeon and the pomeron are quasi-particles which correspond 
to sums of certain Feynman diagrams in the Regge limit (|t| ^ s) |2^. The cross sections for 
reggeon and pomeron exchange have different but universal energy dependences and account 
for all of the total cross section at high energy. There are many different Regge theory- 



based cross section parametrizations possible. Here we use a recent cross section fit |^ 
based on the Donnachie-Landshoff model 



/ 2\ -0-34 / 2\ 0.095 

s — mt^\ s — mt\ 

dreg oc I ^1 CTpom OC I ^\ , (8) 

with the reference scale Sq = 1 GeV^. 

Concerning high energy processes, it is convenient to distinguish between diffractive and 
non-diffractive interactions. Diffractive interactions are characterized by the production of 
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very few secondaries along the direction of the incoming particles. They correspond to the 
quasi-elastic exchange of a reggeon or pomeron between virtual hadronic states the photon 
(mainly the vector mesons c^, and 0) and the nucleon. Because of the spacelike nature of 
the interaction, the angular distribution is strongly forward peaked, and can be parametrized 
by Eq. (|^) with an energy- dependent slope bdis = 6GeV~^ + 0.5GeV~^ ln(s/so) At high 
energies, the cross section of diffractive interactions is approximately a constant fraction of 
the total cross section. The relative contribution of the different vector mesons is predicted 
by theory ||2^ to p° : u; = 9 : 1. The diffractive production of or heavier mesons is neglected 
in SOPHIA. 

Our treatment of non-diffractive multiparticle production is based on the Dual Parton 



Model |]2^. This model can be considered as a phenomenological realization of the expansion 
of QCD for large numbers of colors and flavours [^, ^ in connection with general ideas of 
Duality and Gribov's Regge theory [^, It provides a well developed basic scheme for 



the simulation of high energy hadronic interactions. The model can be visualized as follows: 

(i) the incoming nucleon and photon are split into colored quark and diquark constituents, 

(ii) in the course of the interaction these constituents exchange color quantum numbers, and 

(iii) confinement and the color field forces result in color strings which fragment to hadrons. 
To relate the contributions of reggeon and pomeron exchange to parton configurations, 

we use the correspondence of their respective amplitudes to certain color flow topologies |^ 
which are shown in Figs. |l] and The pomeron exchange topology involves the formation of 
two color neutral strings, while in case of a reggeon topology only one string is stretched from 
the diquark to the quark of the photon. The quark and diquark flavors at the string ends 
are determined by the the spin and valence flavor statistics for the nucleon. For photons the 
charge difference between u and d quarks increases the probability that the photon couples 
to a u-u pair instead of a d-d pair. In the model we use the theoretically predicted ratio of 
4:1 between these two combinations. 

The longitudinal momentum fractions x, 1 — x of the partons connected to the string 
ends are given by Regge asymptotics ||33|, Bl, p5|, BH]. One gets for the valence quark (x) and 



diquark (1 — x) distribution inside the nucleon 

Pix) ^{1 - X^-' (9) 



X 



and for the quark antiquark distribution inside the photon 



p{x) ~ , / . (10) 
ixil — x) 



The relatively small transverse momentum of the partons at the string ends are neglected. 

In the string fragmentation process, the kinetic energy of the initial partons is reduced 
by creating new quark- antiquark pairs which are color field-connected to the parent partons. 
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a) b) 

Figure 1: Color flow picture of (a) a pomeron exchange graph and (b) the final state topology 



q 












qq 






a) b) 

Figure 2: Color flow picture of (a) a reggeon exchange graph and (b) the final state topology 

This process continues until the available kinetic energy drops below the particle production 
threshold causing the newly produced quarks to combine with the valence quarks to form 
hadrons (see, for example, ||37|). 



4 Implementation 



In this section, all numerical expressions for cross sections are measured in units of /ibarn, 
unless noted otherwise. 



4.1 Method of cross section parametrization 

The basic models of photohadronic interaction processes described in Sect. ^ are used to ob- 
tain robust theoretical predictions used for the parametrization of cross sections and final 
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state distribution functions. For theoretically unpredictable parameters we use, if possible, 
the estimates given in the Review of Particle Properties (RPP) ||22|. Remaining parameters 
are determined in fits to available exclusive and inclusive data on 'jp and 7n interactions, 
as compiled in standard reference series ([^, and references therein). Since the parameters 
published in the RPP generally allow some variations within a given error range, these pa- 
rameters are then optimized in comparison to data in an iterative process, until a reasonable 
agreement with the data for a large set of interaction channels is obtained. 



This method has been previously described by Rachen |]I9|, but was considerably im- 
proved in the development of SOPHIA. It provides a minimum bias description of photo- 
production, which reproduces a large set of available data while reducing a possible bias 
due to data selection, since data are used only to fix model parameters. Considering the in- 
tended applications of SOPHIA for (i) astrophysical applications and (ii) the determination 
of background spectra in high energy experiments, we put particular emphasis on a good 
representation of inclusive cross sections and average final state properties in a wide range of 
interaction energies, while a good representation of complex exclusive channels is generally 
not expected. 



4.2 Resonance production 

Using Eq. (Q), the contribution to the cross section from a resonance with mass M and width 
F can be written as a function of the NRF photon energy e' as 

^^') = ^(s-MT + P. • ^^^^ 

The reduced cross section (Tq is entirely determined by the resonance angular momentum 
and the electromagnetic excitation strength b^. We selected all baryon resonances listed in 
the RPP with certain existence (overall status: ****) and a well determined minimal photo- 
excitation strength of > 10~^ for either the prf or the n7 excitation. The resonances 
fulfilling these criteria and their parameters, as implemented in SOPHIA after iterative op- 
timization, are given in Table H. The phase-space reduction close to the Nn threshold is 
heuristically taken into account by multiplying Eq. ( [Tl| ) with the linear quenching function 
Qf(e'; 0.152, 0.17) for the A(1232)-resonance, and with Qf(e'; 0.152, 0.38) for all other reso- 
nances. The function Q{{e'; e[^,w) is defined in Appendix |^. The quenching width w has 
been determined from comparison with the data of the total cross section, and of the ex- 
clusive channels pn^, nir^ and A^+vr" where most of the resonances contribute. The major 
hadronic decay channels of these baryon resonances are A^vr, Att, and N p] for the A^(1535), 
there is also a strong decay into Nr], and the A^(1650) contributes to the KK channel. The 
hadronic decay branching ratios b^ are all well determined for these resonances and given in 
the RPP. However, a difficulty arises from the fact that branching ratios can be expected 
to be energy dependent because of the different masses of the decay products in different 
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Table 2: Baryon resonances and their physical parameters implemented in SOPHIA (see 
text). Superscripts and in the parameters refer to pj and wj excitations, respectively. 
The maximum cross section, (Tmax = 4m^M^(To/(M^ — mj^)"^, is also given for reference. 



resonance 




M 




r 


10^6+ 


4 


a+ 

max 


lO^feO 




max 


A(1232) 


1 


231 





11 


5.6 


31.125 


411.988 


6.1 


33.809 


452.226 


N{UAO) 


1 


440 





35 


0.5 


1.389 


7.124 


0.3 


0.831 


4.292 


A^(1520) 


1 


515 





11 


4.6 


25.567 


103.240 


4.0 


22.170 


90.082 


iV(1535) 


1 


525 





10 


2.5 


6.948 


27.244 


2.5 


6.928 


27.334 


N{1650) 


1 


675 





16 


1.0 


2.779 


7.408 


0.0 


0.000 


0.000 


N{m5) 


1 


675 





15 


0.0 


0.000 


0.000 


0.2 


1.663 


4.457 


N{1680) 


1 


680 





125 


2.1 


17.508 


46.143 


0.0 


0.000 


0.000 


A(1700) 


1 


690 





29 


2.0 


11.116 


28.644 


2.0 


11.085 


28.714 


A(1905) 


1 


895 





35 


0.2 


1.667 


2.869 


0.2 


1.663 


2.875 


A(1950) 


1 


950 





30 


1.0 


11.116 


17.433 


1.0 


11.085 


17.462 



branches. In SOPHIA, we consider all secondary particles, including hadronic resonances, 
as particles of a fixed mass. This implies that, for example, the decay channel Att is ener- 
getically forbidden for ^/s < rriA + '^vr ~ 1-37 GeV. To accommodate this problem, we have 
developed a scheme of energy dependent branching ratios, which change at the thresholds 
for additional decay channels and are constant in between. The requirements are that (i) 
the branching ratio be = for e' < e[-^ j,, and (ii) the average of the branching ratio over 
energy, weighted with the Breit-Wigner function, correspond to the average branching ratio 
given in the RPP for this channel. For all resonances, we considered not more than three 
decay channels leading to a unique solution to this scheme. No fits to data are required. In 
practice, however, the experimental error on many branching ratios allows for some freedom, 
which we have used to generate a scheme that optimizes the agreement with the data on 
different exclusive channels. 

The hadronic branching ratios are given in Tab. ^ in Appendix To obtain the contri- 
bution to a channel with given particle charges, e.g. A~^~^7r~ , the hadronic branching ratio 
^Att has to be multiplied with the iso-branching ratios as given in Tab. 0. We note that with 
the parameters b^, be and 6iso, the resonant contribution to all exclusive decay channels is 
completely determined. 

The angular decay distributions for the resonances follow from Eq. (|^). In SOPHIA, the 
kinematics of the decay channels into Nir is implemented in full detail (see Tab. For other 
decay channels, we assume isotropic decay according to the phase space. Furthermore, there 
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might be some mixing of the different scattering angular distributions since the sampled 
resonance mass, in general, does not coincide with its nominal mass. This effect is neglected 
in our work. Instead, we use the angular distributions applying to resonance decay at its 
nominal mass M. 

Table 3: Angular distribution probability functions for A^vr decay of resonances considered 
in SOPHIA. The resonances N(1535), N(1650) and N(1440) decay isotropically. 



resonance 


^(cos X*) 


A(1232) 


0.636263 - 0.408790 cos2 x* 




Ar0(1520) 


0.673669 - 0.521007 cos2 x* 




Ar+(1520) 


0.739763 - 0.719288 cos2 x* 




A^0(1675) 


0.254005 + 1.427918 cos^ x* - 1-149888 cos^ 


X* 


A^+(1680) 


0.189855 + 2.582610 cos^ x* - 2.753625 cos^ 


X* 


A(1700) 


0.450238 + 0.149285 cos^ x* 




A(1905) 


0.230034 + 1.859396 cos^ x* - 1-749161 cos^ 


X* 


A(1950) 


0.397430 - 1.498240 cos^ x* + 5.880814 cos^ 


X* -4.019252 cos^ X* 



The two decay products of a resonance may also decay subsequently. This decay is 
simulated to occur isotropically according to the available phase space. 

4.3 Direct pion production 

The cross section for direct meson production, unlike those of resonances, is not completely 
determined by well known parameters. The low and high energy constraints suggest the 
phenomenological parametrization 

t^dir(e') = amaxPl(e'; e^^, e'^^^, 2) , (12) 

where the function Pl(e'; e^j^, e'^^^, a) approaches zero for e' = e[^, goes through a maximum 
at e' = e'^^^ and follows an asymptotic behaviour oc (e')^". The definition of this function is 
given in Appendix ^ 

In SOPHIA, we consider explicitly direct channels with charged pion exchange which are 
dominating at low energies. The selection is further constrained by the fact that sufficient 
data are only available for the channels nir^, w-f —>■ pvr", and —>■ A^+vr^. We 

note that proton and neutron induced direct reactions are strictly isospin-symmetric. Both 
proton and neutron data sets (when available) can be used in the fitting procedure. The high 



energy data fits on the Avr and A^vr channel, i.e. cTtt ~ 18(e') ^ and cta ~ 26.4(e') ^ |^ 
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for e' > 1, were primarily used to fix (Xmax, while a best fit of e^^^^ was obtained by comparing 
with the residuals of the low energy data after subtracting the resonance contribution. The 
adopted cross sections are 

a^^(e') = 92.7Pl(e^0.152,0.25,2) + 40expf- ^ Q 002 ) ^^^^ 

^^"""Pl^ 0.002 ) ' 
aA.(e') = 37.7Pl(e';0.4,0.6,2) . (14) 

The two Gaussian-shaped functions included in the direct A^vr cross section have been added 
to improve the representation of the total cross section in the energy region 0.152 GeV < 
e' < 0.4 GeV, where otherwise only the well constrained A(1232) resonance contributes 
significantly. For p^- (wy-) interactions ctatt contributes to the A^+vr" (A"'"7r~) and A^vr"*" 
{A~iT~^) final states with a ratio 3:1 according to isospin combination rules (see Tab. |l]). 

By comparison with the total cross section data we find that the resonant and direct in- 
teraction channels account for all of the total interaction cross section below the 37r threshold 
at e' ~ 0.51 GeV. Above this threshold, and below the threshold for diffractive interactions 
at e' ~ 1 GeV, where high energy processes set in, we find a residual cross section which can 
be parametrized as 

aif = 80.3(60.2)Qf(a;; 0.51; 0.1)(e')^°-^^ , (15) 

where the number 60.2 given in brackets belongs to n7-interactions while the number 80.3 
refers to p7-collisions. The normalization cross section and the quenching width has been 
determined by a minimization method to the total cross section data for p7 (727) inter- 
actions after subtraction of the respective resonant and direct contributions. By analogy, 
the power law index for this contribution is taken from the high energy parametrization for 
reggeon exchange (note that e' oc s — mf^). Physically, this cross section represents the joint 
contribution of all t-channel scattering processes at low energies not considered so far. This 
is in principle similar to interactions at high energies. Consequently, we use an adapted 
string fragmentation model to simulate this contribution, and refer to it as low energy frag- 
mentation hereafter. 

4.4 High energy multipion production 

In SOPHIA, we assume that the cross sections for diffractive and non- diffractive high energy 
interactions are proportional to each other at all energies. This assumption fixes the threshold 
for high energy interactions to the threshold of the Np final state, which is nominally at 
e' ^ 1.1 GeV. Because of the large width of the p there should be some contribution also at 
lower energies. From comparison with exclusive data of the Np final state, and the residuals 
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of the total cross section data, with the sum of the contributions of all low energy channels, 
we find a common threshold for high energy interactions of e^hhigh — 0-85 GeV. 

We restrict the diffractive channel to the non-resonant production of Np and Nu final 
states, for which we assume the theoretically predicted relation cXp = Qcr^^. The ratio between 
diffractive and non-diffractive interactions is derived from the comparison with exclusive Np 
data and total cross section data at high energy, 

CTdifT = 0.15 O-frag . (16) 

For the parametrization of cifrag, we use the power law representations of the reggeon and 
pomeron exchange cross section at high energies, and multiply them by an exponential 
quenching function 1 — exp ( [e(.jj j^jgi^ — e']/a). The relative contributions of the reggeon and 
pomeron cross sections, and the quenching parameter a have been determined by a iterative 
minimization method with respect to the total p7 {nj) cross section data after subtraction 
of all low energy contributions. We find 



Cfragle 



exp 



e' -0.85' 
0.69 



[28.8(26.0) (e')-°-'^ + 58.3{eT ] , (17) 



where we have used the high-energy behaviour given by Eq. (|^). 

The string fragmentation is done by the Lund Monte Carlo Jetset 7.4 ||3^, 0]. This 
program is well suited for string fragmentation at high energies. Since for our purposes also 
strings with rather small invariant masses have to be hadronized, several parameters of this 
fragmentation code had to be tuned to obtain a reasonable description also at low energies. 
Furthermore, in order to avoid double counting, all final states identical to the processes 
already considered by resonance production and direct interactions are rejected (note that 
this is not the case for low-energy fragmentation). 

4.5 Initial state kinematics and photon radiation fields 

The probability for interaction of a proton of energy E^^f with a photon of energy e from a 
radiation field with the photon density n{e) reads 

1 n(e) /"Smax 

where R{En) is the interaction rate as given in Eq. (|^), where also Sth and Smax are defined. 
For a fixed nucleon energy, the CMF energy is sampled from the distribution 

Vis) = <^-\s~m%)aN^{s) (19) 

with $ = Js^'''' ds{s — m\)ai^.y{s). The interaction angle follows from 

^ose = \{'%^ + i\. (20) 



(3 \ 2ENt 
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Currently black body, power law, and broken power law radiation spectra are imple- 
mented in SOPHIA. The photon density n{e) for a blackbody radiation field of temperature 
T is given in natural units by 

^(^) = A T^^-T' (21) 

7r^exp(^) - 1 

where k is the Boltzmann constant. For a power law photon spectrum the photon density is 
given by n{e) — e~'^. The broken power law photon spectrum is given by 

n(e) = e-"i for e < (22) 
n{e) = e^^-"ie-"2 for e > (23) 

where is the break energy, and ai, «2 are the photon indices below and above the break 
energy, respectively. Note that no absolute normalization of n(e) is necessary since it cancels 
in the definition of V{e). 

5 Structure of the program 

5.1 Source code 

The SOPHIA source code consists of several files which contain a number of routines. 

sophia.f main program containing the routines which organize the various tasks to be per- 
formed. Furthermore, the input is handled here and some kinematic transformations 
needed as input to several routines are performed. 

initial. f initialization routine for parameter settings. 

sampling, f collection of routines/functions needed for samphng the CMF energy squared s 
and the photon energy e in the observer frame. 

eventgen.f event generator for photomeson production in jrf and ^7 collisions. 

output. f contains output routines. 

5.2 The event generator 

The simulation of the final state is performed by the photopion production event generator 
EVENTGEN. Together with the initialization routine INITIAL, EVENTGEN can be used separately 
for Monte Carlo event simulation. The user has to give the nucleon code number LO, energy 
EO (in GeV) , the photon of energy e (in GeV) , and the angle 9 in degrees. 

EVENTGEN is structured as follows. First the momenta of the incident particles are Lorentz- 
transformed into the CMF of the interactions. The cross section (in /xbarn) is calculated 
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in the function CROSSECTION. The cross sections for the various channels considered in this 
code determine the distribution for the probabihty of a certain process. The samphng of 
a process (resonance decay, direct channel, diffractive scattering, fragmentation) at a given 
NRF energy of the photon is carried out in the routine DEC_INTER3. 

For the resonance decay we sample the resonance at this energy using the Breit-Wigner 
formula as a probability distribution for a specific resonance (performed in DEC_RES2). Its 
branching ratios define the decay mode (in DEC_PR0C2). The subsequent two-particle decay 
in the CM frame is carried out in RES_DECAY3, and then decays of all unstable particles are 
carried out in the Sybill routine DECS IB. 

Secondary particle production is simulated in GAMMA_H for direct and multiparticle pro- 
duction. 

Finally, all final state particles are Lorentz-transformed back to the LF. 

The output of the final states is organized in the common block /S_PLIST/ P (2000, 5) , 
LLIST(2000) , NP, Ideb. Here the array P(i, j) contains the 4-momenta and rest mass of 
the final state particle i in cartesian coordinates (P(i,l) = P^, P(i,2) = Py, P(i,3) = P^, 
P(i,4) = energy, P(i,5) = rest mass). LLIST() gives the code numbers of all final state 
particles and NP is the number of stable final state particles. 

5.3 Input/Output routines 

Using the standard main program, the user provides the following input parameters. 

• EO = energy of incident proton (in GeV), or 

• Emin, Emax = low/high energy cutoff of an energy grid of incident protons (in GeV) 

• LO = code number of the incident nucleon (LO = 13: proton, LO = 14: neutron) 

• ambient photon field: 

— blackbody spectrum: TBB = temperature (in K) 

— straight /broken power law spectrum: ALPHAl, ALPHA2 = power law indices, 
EPSMIN = low energy cut off (in eV), EPSMAX = high energy cut off (in eV), EPSB 
= break energy (in eV) 

• NTRIAL = number of inelastic interactions 

• NBINS — number of bins for output particle spectra (< 200 bins) 

• DELX = stepsize of output particle spectra 



17 



For the calculation of the incident particle momenta we assume that the relativistic nucleon 
is moving along the positive 2;-axis. 

The output is organized as follows. All the energy distributions (l/Ai'eve)'^A^part/c?logx 
of produced particles are given with logarithmically equal bin sizes in the scaling variable 
X = i?part/EO. Here A'eve denotes the number of simulated inelastic events and A^part is the 
number of secondary particles of a certain kind. The spectra of photons, protons, neutrons, e- 
neutrinos, e-antineutrinos, /i-neutrinos and /i-antineutrinos are considered separately. They 
are stored in a file with name xxxxxx.particle with xxxxxx = input name (chosen by the 
user) : 

particle = 'gamma' — > 7 spectrum 
'muneu' — > spectrum 
'muane' — ^ z/^ spectrum 

'e_neu' — > Ve spectrum 

'eaneu' — > z/g spectrum 

'elect' — > e~ spectrum 

'posit' — > e+ spectrum 

'proto' — > p spectrum 

'neutr' — > n spectrum 
The structure of a typical output file is: 

1. line: high/low energy cutoff of incident nucleon energy grid, 

number of energy bins of incident nucleon spectrum (=NINC), 
TBB or ALPHAl, ALPHA2, EPSMIN, EPSB, EPSMAX, incident particles 

2. line: 1. number: energy of incident nucleon 

2. number: first number (=a) of non-zero energy bin of 

particle spectrum 

3. number: last number (=6) of non-zero energy bin of 

particle spectrum 

3. line: particle spectrum 

between a . . .b 



6 Comparison to data 
6.1 Total cross section 

Fig. 1^ (upper panel) shows the total cross section for 7p-interactions with the various con- 
tributions considered in SOPHIA. For simplicity, we show both fragmentation contributions 
together (low energy fragmentation and non- diffract ive multipion production). The resonant 
process jp —>■ A"*" —>■ ir^p is the only one kinematically possible directly at the particle pro- 
duction threshold. Above the n^n threshold at very low energies (e' < 0.25 GeV), the direct 
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Figure 3: The total cross section (solid line) with the contributions of baryon resonances 
(dotted lines), direct pion production processes (dashed line) non-diffractive multipion pro- 
duction (dash-dotted line), and diffractive scattering (lower dash-triple-dotted line). Data 

Bottom panel: Residuals of the total cross 
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section data to the sum of all partial cross sections implemented in SOPHIA; the line shows 
the average over 10 neighboring data points. 
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channel 7p n^n constitutes the largest contribution. 

To assess the quality of the cross section parametrization, the differences between the 
experimental data on the total 7p cross section and the cross section fit are shown in Fig. ^ 
(lower panel). The total 'jn cross section is overall similar, except at energies of about 
a/s ^ 1.680 GeV (e' ~ 1.035 GeV), where it is considerably smaller due to the different 
excitation strengths of the resonances at this energy. 

6.2 Exclusive cross sections 

Figs. § to compare the output of SOPHIA with the data on specific final states as a function 
of the interaction energy. Such comparisons are important for models that aim to represent 
correctly photohadronic interactions over a wide energy range. 

Fig. ^ compares the total cross sections for vr^p and vr+n production with experimental 
data. The major contributions come from the A(1232) resonance and the direct channel 
together with minor contributions from other resonances. The agreement with data in the 
threshold region is of great importance for many astrophysical applications where this is the 
dominating energy range in the case of steep proton and ambient photon spectra. 

Fig. ^ compares the calculated and measured cross sections for final states involving 
charged and neutral pions. The general agreement of the model results with data is quite 
good, although there are some energy ranges that show minor deviations. It is important to 
note that it is difficult to fit exactly the experimental data without any detailed knowledge 
of the experimental setups and acceptance constraints, especially in cases like 7J9 — > tt'^tt'p 
+ neutrals, where the final state is not well defined. 

The number of inelastic 7p events with 1,3, and 5 charged particles in the final state are 
shown in Fig. ^. This comparison is very sensitive to the description of multipion production 
processes as well as to the smooth transition between different particle production processes. 
It shows that the different channels are reasonably well modeled in SOPHIA. 

Fig. ^ shows reaction cross sections for final states with equal numbers of pions, but 
having a different isospin of the produced nucleon, together with fixed-target data. This 
comparison is important for the correct simulation of the ratio of the proton to neutron 
numbers produced in collisions. Again the model results agree well with data. 

6.3 Inclusive distributions 

The rapidity distribution tests the kinematic of our simulations. The rapidity of a final state 
particle with energy E is defined as 
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Figure 4: Total cross section of 'jp — > vr^p and 7p — > vr+n. Data are from p9| , pO| , [5T| , p2 



53, 54 . 



where py is its momentum component along the direction of the incoming particle. The 
transverse mass ■m± follows from m\ = E"^ — p^y Rapidity is additive under Lorentz trans- 
formations, which keeps the rapidity distribution invariant under such transformations. 
Moffeit et al. have measured the rapidity distribution for the interaction 7p — > 7r^ + 
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anything at three different beam energies e' = 2.8, 4.7 and 9.3 GeV. In Fig. |^ we compare 
the calculated rapidity distribution to data. The agreement in the width and the height of 
the distributions is good. 

6.4 Multiplicities 

In astrophysical environments the production of neutrinos results mainly from the decay 
of charged secondary pions (7r+ — > e'^v^v^Ve^ tt" — > e~z/^z/^z/e). Most of the gamma rays 
are produced via neutral pion decay (vr^ — >■ 77) and synchrotron/Compton emission from 
emanating relativistic e~/e+. Pion multiplicities (see Figs. |10D are therefore instructive to 
understand the energy dissipation of the initial energy among the v- and 7-component. For 
non-astrophysical applications the charged pion multiplicity is a basic interaction parameter 
that presents a cumulative measure of many interaction channels. 

In the resonance region the maximum of the neutral pion multiplicity is reached at the 
A(1232)-resonance (see Fig. At threshold neutral pion production is strongly suppressed 
in favour of vr^-production (see Fig. ITUp due to the dominance of the direct channel. Mul- 
tiplicities are approximately growing as ~ s^/'^ in the multipion production region at low 
energies. The obtained multiplicity distributions from our simulations are in agreement with 
the data from Moffeit et al. |0 . 



By counting the number of protons and neutrons produced in 7p interactions, one can 
define a proton-to-neutron ratio. The SOPHIA prediction on the energy dependence of this 
ratio is shown in Fig. The p/n ratio reaches ~ 2.2 at high energy which can be contrasted 
to an experimentally estimated value of about 3.8 found by Meyer The ratio derived in 



|68| is essentially based on the same data as those our model is compared to. Meyer estimated 
the experimentally unknown cross sections by isospin symmetry arguments. In our case these 
cross sections are predicted by the Monte Carlo simulation. We conclude from this that the 
unmeasured cross sections are important for the p/n ratio, and that the difference between 
the values 2.2 and 3.8 reflect the uncertainty due to the limited experimental data available. 



7 Conclusions 

A newly developed Monte Carlo event generator SOPHIA has been presented. It simulates 
the interactions of nucleons with photons over a wide range in energy. The simulation of the 
final state includes all interaction processes which are relevant to astrophysical applications. 
SOPHIA contains tools, such as for sampling the photon energy from different ambient soft 
photon spectra, and the nucleon-photon interaction angle, that are needed for such appli- 
cations. As an event generator, SOPHIA uses all available information on the interaction 
cross section, the final state particle composition, and kinematics of the interaction processes 
as provided by particle physics. Comparison with the available accelerator data shows that 
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SOPHIA provides a good description of our current knowledge of photon-nucleon interac- 
tions. 
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Figure 8: Rapidity distribution for 7p — > tt +anything at beam energies e' = 2.8GeV 
(dotted line), 4.7 GeV (solid line) and 9.3 GeV (thick solid line). Data are from Moffeit et 
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Figure 10: Left figure: vr^-multiplicity for 7p-interactions. Right figure: Charged particle 
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Figure 11: Proton-to-neutron ratio for 7p-interactions as simulated by SOPHIA. At high 
energy this ratio is ~ 2.2. 
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A Definition of functions 

The functions Pl(e', e[^, e^g^, a) and Qf (e', e[^, w) are defined by 

PKe', 4, 4a., «) = ( /~_';^ (25) 

max th max 

for e' > e[^, and Pl(e', e[^, e'^^, a) = otherwise, A = ae'^^/e[^, and 

Qi(e',e[^,w)^(e'-e'J/w (26) 
for e[^<e' <w + e^^, Qf (e', e^^, w) = for e' < e^^, and Qf (e', e^^, w) = 1 for e' < w + e^^- 
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B Resonance branching ratios 



Table 4: Decay channels and branching ratios for resonance decays, following the scheme 
described in the text. The second row of each table gives the fractional contribution of 
the decay probability for several given energy ranges. The average branching ratio for each 
channel is given in the last column of each table. They are generally consistent with the 
average values quoted in the RPP p2|, except for the Np decay of the iV(1520), where 
bNp < 6% due to the low fractional contribution of this energy range whereas 15%— 25% are 
quoted in the RPP. 



A(1232) 


all X 


total 


fraction 


100% 




Nvr 


100% 


100% 



N(1440) 


x<0.54 x>0.54 


total 




fraction 


34% 66% 








Nvr 


100% 50% 


67% 






Avr 


50% 


33% 




N(1520) 


x<0.54 0.54<x<1.09 x>1.09 


total 


fraction 


9% 


85% 


6% 




Nvr 


100% 


50% 


0% 


52% 


Att 




50% 


0% 


42% 


Np 






100% 


6% 


N(1535) 


x<0.71 x>0.71 


total 




fraction 


31% 69% 








Nvr 


100% 25% 


45% 






N?7 


75% 


55% 




N(1650) 


x<0.54 0.54<x<0.91 x>0.91 


total 


fraction 


5% 


26% 


69% 




Nvr 


100% 


85% 


70% 


75% 


Avr 




15% 


15% 


14% 


AK 






15% 


11% 
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N(1675) 


x<0.54 


x>0.54 


total 


fraction 


5% 


95% 




Nvr 


100% 


40% 


42% 


Att 




60% 


57% 



N(1680) 


x<0.54 


0.54<x<1.09 


x>1.09 


total 


fraction 


4% 


64% 


32% 




Ntt 


100% 


65% 


55% 


64% 


Att 




35% 


0% 


22% 


Np 






45% 


14% 


A(1700) 


x<0.54 


0.54<x<1.09 


x>1.09 


total 


fraction 


9% 


52% 


39% 




Ntt 


100% 


0% 


0% 


14% 


Att 


- 


100% 


20% 


55% 


N/9 






80% 


31% 


A(1905) 


x<1.09 


0.54<x<1.09 


x>1.09 


total 


fraclicni 




21% 


73% 




Ntt 


100% 


40% 


0% 


14% 


Att 




60% 


0% 


13% 


N/9 






100% 


73% 


A(1950) 


x<0.54 


0.54<x<1.09 


x>1.09 


total 


fraction 


4% 


15% 


81% 




Ntt 


100% 


60% 


30% 


37% 


Att 




40% 


40% 


39% 


N/9 






30% 


24% 



35 



C Compilation of routines/functions 



function BREITWIGNER(SIGMA_0, GAMMA, DMM, EPS_PRIME): 

calculates cross section (in /xbarn) of a resonance with width GAMMA (in GeV) , mass DMM 
(in GeV), maximum cross section SIGMA_0 (in //barn) and NRF energy of the photon 
(in GeV) according to the Breit-Wigner formula 

subroutine CROSSDIR(EPS_PRIME): 

collection of functions (SINGLEBACK, TWOBACK) which calculates the cross section (at 
the NRF energy EPS_PRIME of the incident photon) of the direct channel (not isospin- 
corrected) 

function CROSSECTION(EPS_PRIME,NDIR): 

computes cross section (in //barn) for A'^7-interaction at a given energy EPS_PRIME 
(=photon energy in GeV in proton's rest frame); depending on the control variable 
NDIR it returns the total cross section (NDIR=3) or only a certain part of the cross 
section (NDIR=1: total resonance cross section, NDIR=4: direct channel cross section 
NDIR=5: fragmentation cross section, NDIR=11-19: individual resonance cross sections) 

subroutine DEC_INTER3 (EPS_PRIME , IMODE) : 

returns reaction mode: decay of resonance (IM0DE=6), direct pion production (IM0DE=2 
for Nt: final states, IM0DE=3 for Att final states), fragmentation in resonance re- 
gion (IM0DE=5), diffractive scattering (IM0DE=1 for Np final states, IM0DE=4 for Nuj 
final states) and multipion production/fragmentation (IM0DE=0) at a given energy 
EPS_PRIME (in GeV) 

subroutine DEC_PR0C2 (EPS_PRIME , IPROC , IRANGE , IRES , LO) : 

returns in IPROC the decay mode for a given resonance IRES at energy EPS_PRIME (in 
GeV) for incident nucleon with code number LO; IRANGE is the number of energy inter- 
vals corresponding to the (energy dependent) branching ratios of a specific resonance 

subroutine DEC_RES2 (EPS_PRIME , IRES , IRESMAX , LO) : 

returns sampled resonance number IRES (IRES = 1 . . . IRESMAX) at energy EPS_PRIME 
(in GeV) for incident nucleon with code number LO 

subroutine DECSIB: 

decay of unstable particles 

function FUNCTS(S): 

calculates distribution of the squared CMF energy S (in GeV^) 

subroutine GAMMA_H(EO,LO, IMODE): 

interface routine for hadron-7 collisions in CM frame; 
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EO is CMF energy of the A'"7-system; first particle (7 or hadron N) goes in +2;-direction; 
final state consists of protons, neutrons, gamma, leptons, neutrinos 

routine INITIAL (LO): 

initialization routine for parameter settings; 
to be called at before calling EVENTGEN 

routine LUND_FRAG(SQS,NP): 

interface to Jetset fragmentation of system with CMF energy SQS (in GeV); 
NP = number of secondaries produced 

function PHOTD (EPS , TBB) : 

returns photon density (in photons/cm^/eV) at energy EPS (in eV) of blackbody radi- 
ation of temperature TBB (in K) 

subroutine PROBANGLE ( IRES , LO , ANGLESCAT) : 

probability distribution for scattering angle of given resonance IRES and incident nu- 
cleon with code number LO ; 

ANGLESCAT is the cosine of scattering angle in CMF frame 

function PROB_EPSKT(EPS) : 

calculates probability distribution for thermal photon field with temperature TBB (in 
K) at energy EPS (in eV) 

function PROB_EPSPL(EPS) : 

calculates probability distribution at energy EPS (in eV) for power law radiation n ~ 
£pg-ALPHA between the limits EPSMl . . .EPSM2 (in eV) 

subroutine PROC_TWOPART (LA , LB , AMD , LRES , PRES , COSTHETA , NBAD) : carries out 2-particle 
decay in CM frame of system with mass AMD (in GeV) into particles with code numbers 
LA, LB (stored in array LRESO) wheras COSTHETA is the cosine of the N^-CM frame 
scattering angle ; returns 5-momenta PRESO of the two particles; 
NBAD=1 for kinematically not possible decays, otherwise NBAD=0 

subroutine RES_DECAY3 (IRES , IPROC , IRANGE , S , LO , NBAD) : 

determines decay products for the decay of a given resonance IRES by a given decay 
process IPROC for a incident nucleon with code number LO, and carries out two-particle 
decay of the resonance with squared CMF energy S; code numbers and (CMF) 5- 
momenta of produced particles are stored in arrays LLIST and P, respectively, in a 
common block; 

NBAD=1 for kinematically not possible decays; otherwise NBAD=0 
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subroutine SAMPLE_EPS (EPS , EPSMIN , EPSMAX) : 

samples incident photon energy; if the blackbody temperature is set to TBB=0 the 
photon energy is sampled from a power law distribution with index ALPHA between given 
photon energies EPSMIN, EPSMAX (in eV); otherwise it is sampled from a distribution 
for a blackbody radiation with temperature TBB (in K) 

subroutine SAMPLE_S (S ,EPS) : 

samples the total CM frame energy S (in GeV^) for a given photon with energy EPS 
(in GeV) and proton with energy EO (in GeV) 

subroutine SCATANGLE (ANGLESCAT , IRES , IPROC , INC) : 

samples the cosine of the scattering angle ANGLESCAT for a given resonance IRES and 
incident nucleon INC 
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